library(edgeR)
library(limma)
library(SingleCellExperiment)
library(SingleR)
library(Seurat)
library(biomaRt)
library(dplyr)
library(reshape2)
library(patchwork)
library(ComplexHeatmap)
library(ggplotify)
library(gridExtra)
library(grid)
library(viridis)
Loading required package: limma
Loading required package: SummarizedExperiment
Loading required package: MatrixGenerics
Loading required package: matrixStats
Attaching package: ‘MatrixGenerics’
The following objects are masked from ‘package:matrixStats’:
colAlls, colAnyNAs, colAnys, colAvgsPerRowSet, colCollapse,
colCounts, colCummaxs, colCummins, colCumprods, colCumsums,
colDiffs, colIQRDiffs, colIQRs, colLogSumExps, colMadDiffs,
colMads, colMaxs, colMeans2, colMedians, colMins, colOrderStats,
colProds, colQuantiles, colRanges, colRanks, colSdDiffs, colSds,
colSums2, colTabulates, colVarDiffs, colVars, colWeightedMads,
colWeightedMeans, colWeightedMedians, colWeightedSds,
colWeightedVars, rowAlls, rowAnyNAs, rowAnys, rowAvgsPerColSet,
rowCollapse, rowCounts, rowCummaxs, rowCummins, rowCumprods,
rowCumsums, rowDiffs, rowIQRDiffs, rowIQRs, rowLogSumExps,
rowMadDiffs, rowMads, rowMaxs, rowMeans2, rowMedians, rowMins,
rowOrderStats, rowProds, rowQuantiles, rowRanges, rowRanks,
rowSdDiffs, rowSds, rowSums2, rowTabulates, rowVarDiffs, rowVars,
rowWeightedMads, rowWeightedMeans, rowWeightedMedians,
rowWeightedSds, rowWeightedVars
Loading required package: GenomicRanges
Loading required package: stats4
Loading required package: BiocGenerics
Loading required package: parallel
Attaching package: ‘BiocGenerics’
The following objects are masked from ‘package:parallel’:
clusterApply, clusterApplyLB, clusterCall, clusterEvalQ,
clusterExport, clusterMap, parApply, parCapply, parLapply,
parLapplyLB, parRapply, parSapply, parSapplyLB
The following object is masked from ‘package:limma’:
plotMA
The following objects are masked from ‘package:stats’:
IQR, mad, sd, var, xtabs
The following objects are masked from ‘package:base’:
anyDuplicated, append, as.data.frame, basename, cbind, colnames,
dirname, do.call, duplicated, eval, evalq, Filter, Find, get, grep,
grepl, intersect, is.unsorted, lapply, Map, mapply, match, mget,
order, paste, pmax, pmax.int, pmin, pmin.int, Position, rank,
rbind, Reduce, rownames, sapply, setdiff, sort, table, tapply,
union, unique, unsplit, which.max, which.min
Loading required package: S4Vectors
Attaching package: ‘S4Vectors’
The following objects are masked from ‘package:base’:
expand.grid, I, unname
Loading required package: IRanges
Loading required package: GenomeInfoDb
Loading required package: Biobase
Welcome to Bioconductor
Vignettes contain introductory material; view with
'browseVignettes()'. To cite Bioconductor, see
'citation("Biobase")', and for packages 'citation("pkgname")'.
Attaching package: ‘Biobase’
The following object is masked from ‘package:MatrixGenerics’:
rowMedians
The following objects are masked from ‘package:matrixStats’:
anyMissing, rowMedians
Attaching package: ‘SingleCellExperiment’
The following object is masked from ‘package:edgeR’:
cpm
Attaching SeuratObject
Attaching package: ‘Seurat’
The following object is masked from ‘package:SummarizedExperiment’:
Assays
Attaching package: ‘dplyr’
The following object is masked from ‘package:biomaRt’:
select
The following object is masked from ‘package:Biobase’:
combine
The following objects are masked from ‘package:GenomicRanges’:
intersect, setdiff, union
The following object is masked from ‘package:GenomeInfoDb’:
intersect
The following objects are masked from ‘package:IRanges’:
collapse, desc, intersect, setdiff, slice, union
The following objects are masked from ‘package:S4Vectors’:
first, intersect, rename, setdiff, setequal, union
The following objects are masked from ‘package:BiocGenerics’:
combine, intersect, setdiff, union
The following object is masked from ‘package:matrixStats’:
count
The following objects are masked from ‘package:stats’:
filter, lag
The following objects are masked from ‘package:base’:
intersect, setdiff, setequal, union
Loading required package: grid
========================================
ComplexHeatmap version 2.8.0
Bioconductor page: http://bioconductor.org/packages/ComplexHeatmap/
Github page: https://github.com/jokergoo/ComplexHeatmap
Documentation: http://jokergoo.github.io/ComplexHeatmap-reference
If you use it in published research, please cite:
Gu, Z. Complex heatmaps reveal patterns and correlations in multidimensional
genomic data. Bioinformatics 2016.
The new InteractiveComplexHeatmap package can directly export static
complex heatmaps into an interactive Shiny app with zero effort. Have a try!
This message can be suppressed by:
suppressPackageStartupMessages(library(ComplexHeatmap))
========================================
Attaching package: ‘gridExtra’
The following object is masked from ‘package:dplyr’:
combine
The following object is masked from ‘package:Biobase’:
combine
The following object is masked from ‘package:BiocGenerics’:
combine
Loading required package: viridisLite
# Set working directory to project root
setwd("/research/peer/fdeckert/FD20200109SPLENO")
# Source files
source("plotting_global.R")
Loading required package: RColorBrewer
# Seurat files
cnt_file <- "data/haemosphere/Haemopedia-Mouse-RNASeq_raw.txt"
meta_file <- "data/haemosphere/Haemopedia-Mouse-RNASeq_samples.txt"
so_file <- "data/object/seurat.rds"
filter_min_cpm = 0.5
filter_min_expressed_samples = 2
normalization_method = "TMM"
adj_p_cutoff = 0.05
min_rows = 50
# Plotting Theme
ggplot2::theme_set(theme_global_set()) # From project global source()
# Settings
options(digits = 4)
Loading required package: ggplot2
# Takes a prcomp pca object and group labels as input to compute a PCA scattern plot of the first two PCA
pca_scattern <- function(pca, group, group_color) {
pca_data <- cbind(pca$x, data.frame(group = group))
pca_scattern <- ggplot(pca_data, aes(x = PC1, y = PC2, fill = group)) +
geom_point(color = "black", pch = 21, size = 3) +
ggtitle("log2(CPM) - z-score") +
xlab(paste("PC1 (", round(100 * (pca$sdev^2)[1] / sum(pca$sdev^2), digits = 1), " %)")) +
ylab(paste("PC2 (", round(100 * (pca$sdev^2)[2] / sum(pca$sdev^2), digits = 1), " %)")) +
scale_fill_manual(values = group_color) +
guides(fill = guide_legend(title = "Group")) +
theme(aspect.ratio = 1,
panel.border = element_blank(),
axis.line = element_line())
return(pca_scattern)
}
top_hm <- function(deg, cnt = NULL, meta = NULL, level = "", contrast = "", top = 50, min_contrast_log2cpm = 2) {
# Filter and sort cnt by meta samples
# Columns of cnt must be corespondend to meta sample id to compute mean cpm by group
cnt <- cnt[, colnames(cnt) %in% meta$sample_id]
cnt <- cnt[, match(meta$sample_id, colnames(cnt))]
# Compute cpm and log2(cpm)
cpm <- edgeR::cpm(cnt, normalized.lib.sizes = FALSE, log = FALSE, prior.count = 0)
log_cpm <- edgeR::cpm(cnt, normalized.lib.sizes = FALSE, log = TRUE, prior.count = 1)
# Define grouping levels to compute the group mean log2 cpm
if(level == "main_labels") {group <- meta$main_labels}
if(level == "fine_labels") {group <- meta$fine_labels}
# Mean group cpm and log2(cpm)
colnames(cpm) <- group
cpm <- lapply(unique(group), function(i) data.frame(X = rowMeans(cpm[, grepl(i, colnames(cpm), fixed = TRUE)])))
cpm <- do.call("cbind", cpm)
colnames(cpm) <- unique(group)
colnames(log_cpm) <- group
log_cpm <- lapply(unique(group), function(i) data.frame(X = rowMeans(log_cpm[, grepl(i, colnames(log_cpm), fixed = TRUE)])))
log_cpm <- do.call("cbind", log_cpm)
colnames(log_cpm) <- unique(group)
# Compute log2fc between contrast and maximum group
max_log2fc <- data.frame(max_log2fc = log2(cpm[, contrast] / apply(cpm[, colnames(cpm) != contrast], 1, max)))
rownames(max_log2fc) <- rownames(cpm)
# filter log_cpm and max_log2fc by deg top hits
log_cpm <- log_cpm[rownames(log_cpm) %in% rownames(deg), ]
max_log2fc <- max_log2fc[rownames(max_log2fc) %in% rownames(deg), , drop = FALSE]
# Sort log_cpm by max_log2fc
log_cpm <- log_cpm[arrange(max_log2fc, desc(max_log2fc)) %>% rownames(), ]
# Filter for min expression in contrast group
log_cpm <- log_cpm[log_cpm[, contrast] > min_contrast_log2cpm, ]
# Filter log_cpm top hits
if (nrow(log_cpm) > top) {log_cpm <- log_cpm[1:top, ]}
# ComplexHeatmap color ramp
color_range <- max(abs(log_cpm))
color_break <- seq(0, color_range, 0.01)
color_ramp <- viridis(length(color_break), option = "magma")
# Heat map
hm <- grid.grabExpr(draw(ComplexHeatmap::pheatmap(
mat = as.matrix(log_cpm),
main = contrast,
fontsize_row = 10,
scale = "none",
cluster_rows = TRUE,
cluster_cols = FALSE,
cellwidth = 12,
cellheight = 12,
clustering_distance_rows = "euclidean",
clustering_distance_cols = "euclidean",
clustering_method = "complete",
show_row_dend = FALSE,
show_rownames = TRUE,
show_colnames = TRUE,
color = color_ramp,
breaks = color_break,
border_color = NA)))
return(hm)
}
# Contrast matrix for one vs average of all other explanatory variables
contrast_one_vs_average <- function(design) {
exp_variables <- colnames(design)
contrast <- list()
for(i in seq_along(exp_variables)) {
design_string <- paste(exp_variables[i], "-", "(", paste(exp_variables[-i], collapse = "+"), ")/", length(exp_variables[-i]), sep = "")
contrast[[i]] <- makeContrasts(contrasts = design_string, levels = exp_variables)
}
return(contrast)
}
# Compute DEG matrix needs a contrast, design and count matrix object
deg <- function(
cnt, design, contrast,
filter_min_cpm = 0.5,
filter_min_expressed_samples = 2,
normalization_method = "TMM",
adj_p_cutoff = 0.05
) {
# Filter and sort cnt by design matrix
cnt <- cnt[, colnames(cnt) %in% rownames(design)]
cnt <- cnt[, match(rownames(design), colnames(cnt))]
# Create a DGEList object from cnt and normalize
cnt <- DGEList(cnt)
cnt <- calcNormFactors(cnt, method = normalization_method)
# Keep only rows where filter_min_expressed_samples or more samples out of group1 and group2 must be greater than filter_min_cpm
cnt <- cnt[rowSums(edgeR::cpm(cnt, normalized.lib.sizes = TRUE, log = FALSE, prior.count = 0) > filter_min_cpm) >= filter_min_expressed_samples, ] # cpm default is normalized.lib.sizes = TRUE and log = FALSE
# Limma fit to fit explanatory variable
fit <- lmFit(voom(cnt, design = design, normalize.method = "none"), design)
# For a linear model fit, compute moderated t-statistics, moderated F-statistic, and log-odds of differential expression by empirical Bayes moderation of the standard errors towards a common value.
eb <- eBayes(contrasts.fit(fit, contrast))
# Fetch topTable
m = topTable(eb, adjust = "fdr", number = Inf) # This should fetch every gene/probe sorted by adj p
m = m[m$adj.P.Val < adj_p_cutoff, ]
# select column subset and remove dots in adj.P.Val colname
m = m[c("logFC", "adj.P.Val","AveExpr")] # missing logFC
colnames(m)[2] = "adjPValue"
return(m)
}
cnt <- read.delim(cnt_file, row.names = 1)
meta <- read.delim(meta_file)
# Select sample by cell_name
cell_name <- c(
"Fob", "MZB", "GCB", "SplPlsC", #B-cells SPLEEN
"BasoBM", #Basophil BM
"cDC2", "pDC", #DC conventional and plasmacytoid SPLEEN
"EoP", "Eo", #Eosinophil progenitor and adult BM
"MEP", "PreCFUE", "CFUE", "EryBlPB", "EryBlPO", "Retic", #Erythrocyte BM
"MonoBM", "Mac", #Monocyte BM and macrophage PERITONEAL
"Mast", #Mast cell PERITONEAL
"MegTPO", #Megakaryocyte cultured from BM
"STHSC", "LSK", "MPP", #Multi potential progenitor BM
"NeutBM", #Neutrophil BM
"NK", #Natural killer cell SPLEEN
"CMP", #Common myeloid Progenitor BM
"GMP", #Granulocyte macrophage progenitor BM
"CLP", #Common lymphoid progenitor BM
"CD4T", "CD8T" #T-cells total SPLEEN and LYMPHNODES
)
meta <- meta[meta$cell_name %in% cell_name, ]
# Filter and order count matrix columns by meta data samples
cnt <- cnt[, colnames(cnt) %in% meta$sample_id]
cnt <- cnt[, match(meta$sample_id, colnames(cnt))]
meta$sample_id <- make.names(paste0(meta$fine_labels, ".", seq(1:nrow(meta))))
colnames(cnt) <- meta$sample_id
# Get ensembl to symbol relationship from biomaRt
mart <- useDataset("mmusculus_gene_ensembl", useMart("ensembl"))
ens_to_sym <- getBM(filters = "ensembl_gene_id", attributes = c("ensembl_gene_id", "mgi_symbol"), values = rownames(cnt), mart = mart)
# Remove duplicates
ens_to_sym <- ens_to_sym[!ens_to_sym$ensembl_gene_id %>% duplicated, ]
ens_to_sym <- ens_to_sym[!ens_to_sym$mgi_symbol %>% duplicated, ]
# Filter data by ensembl id in reference
cnt <- cnt[rownames(cnt) %in% ens_to_sym$ensembl_gene_id, ]
# order ens_to_sym ensembl id by rownames of data
ens_to_sym <- ens_to_sym[match(ens_to_sym$ensembl_gene_id, rownames(cnt)), ]
# Set mgi_symbol as data rownames
rownames(cnt) <- ens_to_sym$mgi_symbol
# Keep only genes with counts greater 0
cnt <- cnt[rowSums(cnt) > 0, ]
# Set factor level for labels
meta$main_labels <- factor(meta$main_labels, levels = names(color$main_labels_haemosphere))
meta$fine_labels <- factor(meta$fine_labels, levels = names(color$fine_labels_haemosphere))
options(repr.plot.width = 10, repr.plot.height = 5)
# create a DGEList object from cnt and normalize
cnt_pca <- cnt[rowSums(edgeR::cpm(cnt, normalized.lib.sizes = FALSE, log = FALSE) > 0.5) >= 2, ]
pca <- prcomp(t(edgeR::cpm(cnt_pca, normalized.lib.sizes = FALSE, log = TRUE)), scale = TRUE, center = TRUE)
pca_1 <- pca_scattern(pca, meta$main_labels, color$main_labels_haemosphere)
pca_2 <- pca_scattern(pca, meta$fine_labels, color$fine_labels_haemosphere)
pca_1 + pca_2
Explanatory variable (factor): Cell lineage (design_cl) or cell type (design_ct)
Indipendent meassurements (character): Cell type with batch information
Statistical model: Means model to get all possible combinations of levels
Design matrix: Coded as model.matrix(~ 0+explanatory_variable) which is here equivalent to model.matrix(~ explanatory_variable)
meta_ml <- meta
meta_ml$main_labels <- make.names(meta_ml$main_labels)
# Make design for means model parameterization
design_ml = model.matrix(~0+as.factor(meta_ml$main_labels))
colnames(design_ml) = unique(meta_ml$main_labels)
rownames(design_ml) = meta_ml$sample_id
# Create contrast
contrast_ml <- contrast_one_vs_average(design_ml)
# Compute deg for all contrast combination
deg_ml <- list()
for(i in seq_along(contrast_ml)) {deg_ml[[i]] <- deg(cnt = cnt, design = design_ml, contrast = contrast_ml[[i]])}
options(repr.plot.width = 20, repr.plot.height = 40)
# Plot heat map of marker genes
hm_ml <- lapply(seq_along(unique(meta_ml$main_labels)), function(i) top_hm(deg = deg_ml[[i]], cnt = cnt, meta = meta_ml, level = "main_labels", contrast = unique(meta_ml$main_labels)[i], top = 50, min_contrast_log2cpm = 2))
gridExtra::arrangeGrob(grobs = hm_ml, ncol = 4) %>% grid::grid.draw()
meta_fl <- meta
meta_fl$fine_labels <- make.names(meta_fl$fine_labels)
# Make design for means model parameterization
design_fl <- model.matrix(~0+as.factor(meta_fl$cell_name))
colnames(design_fl) <- unique(meta_fl$cell_name)
rownames(design_fl) <- meta_fl$sample_id
# Create contrast
contrast_fl <- contrast_one_vs_average(design_fl)
# Compute deg for all contrast combination
deg_fl <- list()
for(i in seq_along(contrast_fl)) {deg_fl[[i]] <- deg(cnt = cnt, design = design_fl, contrast = contrast_fl[[i]])}
options(repr.plot.width = 20, repr.plot.height = 90)
# Plot heat map of marker genes
hm_fl <- lapply(seq_along(unique(meta_fl$fine_labels)), function(i) top_hm(deg = deg_fl[[i]], cnt = cnt, meta = meta_fl, level = "fine_labels", contrast = unique(meta_fl$fine_labels)[i], top = 50, min_contrast_log2cpm = 1))
gridExtra::arrangeGrob(grobs = hm_fl, ncol = 3) %>% grid::grid.draw()
meta_ery <- meta[meta$main_labels == "Ery", ]
meta_ery$fine_labels <- make.names(meta_ery$fine_labels)
# Make design for means model parameterization
design_ery <- model.matrix(~0+as.factor(meta_ery$cell_name))
colnames(design_ery) <- unique(meta_ery$cell_name)
rownames(design_ery) <- meta_ery$sample_id
# Create contrast
contrast_ery <- contrast_one_vs_average(design_ery)
# Compute deg for all contrast combination
deg_ery <- list()
for(i in seq_along(contrast_ery)) {deg_ery[[i]] <- deg(cnt = cnt, design = design_ery, contrast = contrast_ery[[i]])}
options(repr.plot.width = 20, repr.plot.height = 10)
# Plot heat map of marker genes
hm_ery <- lapply(seq_along(unique(meta_ery$fine_labels)), function(i) top_hm(deg = deg_ery[[i]], cnt = cnt, meta = meta_ery, level = "fine_labels", contrast = unique(meta_ery$fine_labels)[i], top = 50, min_contrast_log2cpm = 1))
gridExtra::arrangeGrob(grobs = hm_ery, ncol = 6) %>% grid::grid.draw()
meta_mpp <- meta[meta$main_labels == "MPP", ]
meta_mpp$fine_labels <- make.names(meta_mpp$fine_labels)
# Make design for means model parameterization
design_mpp <- model.matrix(~0+as.factor(meta_mpp$cell_name))
colnames(design_mpp) <- unique(meta_mpp$cell_name)
rownames(design_mpp) <- meta_mpp$sample_id
# Create contrast
contrast_mpp <- contrast_one_vs_average(design_mpp)
# Compute deg for all contrast combination
deg_mpp <- list()
for(i in seq_along(contrast_mpp)) {deg_mpp[[i]] <- deg(cnt = cnt, design = design_mpp, contrast = contrast_mpp[[i]])}
options(repr.plot.width = 20, repr.plot.height = 10)
# Plot heat map of marker genes
hm_mpp <- lapply(seq_along(unique(meta_mpp$fine_labels)), function(i) top_hm(deg = deg_mpp[[i]], cnt = cnt, meta = meta_mpp, level = "fine_labels", contrast = unique(meta_mpp$fine_labels)[i], top = 50, min_contrast_log2cpm = 1))
gridExtra::arrangeGrob(grobs = hm_mpp, ncol = 6) %>% grid::grid.draw()
meta_rpp <- meta[meta$main_labels == "RPP", ]
meta_rpp$fine_labels <- make.names(meta_rpp$fine_labels)
# Make design for means model parameterization
design_rpp <- model.matrix(~0+as.factor(meta_rpp$cell_name))
colnames(design_rpp) <- unique(meta_rpp$cell_name)
rownames(design_rpp) <- meta_rpp$sample_id
# Create contrast
contrast_rpp <- contrast_one_vs_average(design_rpp)
# Compute deg for all contrast combination
deg_rpp <- list()
for(i in seq_along(contrast_rpp)) {deg_rpp[[i]] <- deg(cnt = cnt, design = design_rpp, contrast = contrast_rpp[[i]])}
options(repr.plot.width = 20, repr.plot.height = 10)
# Plot heat map of marker genes
hm_rpp <- lapply(seq_along(unique(meta_rpp$fine_labels)), function(i) top_hm(deg = deg_rpp[[i]], cnt = cnt, meta = meta_rpp, level = "fine_labels", contrast = unique(meta_rpp$fine_labels)[i], top = 50, min_contrast_log2cpm = 1))
gridExtra::arrangeGrob(grobs = hm_rpp, ncol = 6) %>% grid::grid.draw()
# Create SingleR object from Haemosphere RNA-seq data
ref <- SummarizedExperiment(list(counts = cnt, logcounts = edgeR::cpm(cnt, log = TRUE, prior.count = 1)))
ref$label.main <- as.character(meta$main_labels)
ref$label.fine <- as.character(meta$fine_labels)
# Import Seurat test data
so <- readRDS(so_file)
# Seurat to SingleCellExperiment
sce <- SingleCellExperiment(list(counts = so@assays$RNA@counts))
# Predict labels
label_main <- SingleR::SingleR(test = sce, ref = ref, labels = ref$label.main, assay.type.test = "counts", de.method = "classic")
label_fine <- SingleR::SingleR(test = sce, ref = ref, labels = ref$label.fine, assay.type.test = "counts", de.method = "classic")
# Add labels to Seurat object
label_main_meta <- as.data.frame(label_main) %>% dplyr::select(pruned.labels, tuning.scores.first) %>% dplyr::rename(main_labels = pruned.labels, main_delta_score = tuning.scores.first)
label_fine_meta <- as.data.frame(label_fine) %>% dplyr::select(pruned.labels, tuning.scores.first) %>% dplyr::rename(fine_labels = pruned.labels, fine_delta_score = tuning.scores.first)
so <- AddMetaData(so, cbind(label_main_meta, label_fine_meta))
# Set factor level for labels
so$main_labels <- factor(so$main_labels, levels = names(color$main_labels_haemosphere))
so$fine_labels <- factor(so$fine_labels, levels = names(color$fine_labels_haemosphere))
# write.csv(cbind(label_main_meta, label_fine_meta), "data/object/components/annotation/singler_haemosphre_counts.csv")
so <- AddMetaData(so, read.csv("data/object/components/annotation/singler_haemosphre_counts.csv", row.names = 1))
The main labels are not predicted by SingleR! We just use the top annoation for grouping e.g. the top lable for a particular cell might differ.
options(repr.plot.width = 15, repr.plot.height = 8)
hm <- function(data, fill) {
set <- data
# Mean delta score for fine labels per tissue group
data <- dplyr::group_by(data, fine_labels, tissue) %>% dplyr::summarise(mean_delta_score = mean(fine_delta_score))
# Add main_label information for fine_labesl
data <- dplyr::left_join(meta[, c("fine_labels", "main_labels")] %>% unique(), data, by = "fine_labels")
# Reshape to get scores by tissue
data <- reshape2::dcast(data, fine_labels+main_labels~tissue, value.var = "mean_delta_score")
# Set NA scores to 0 for plotting
data$Myeloid[is.na(data$Myeloid)] <- 0
data$Progenitor[is.na(data$Progenitor)] <- 0
data$Myeloid[is.na(data$fine_labels)] <- 0
data$Progenitor[is.na(data$fine_labels)] <- 0
# Set factor levels
data$main_labels <- factor(data$main_labels, levels = names(color$main_labels_haemosphere))
data$fine_labels <- factor(data$fine_labels, levels = names(color$fine_labels_haemosphere))
# Get matrix for heatmap
m <- t(as.matrix(data[, c("Myeloid", "Progenitor")]))
colnames(m) <- data$fine_labels
# Column annotation (not used)
col_anno <- data.frame(data$main_labels)
colnames(col_anno) <- c('Cell lineage')
col_colors <- list("Cell lineage" = color$main_labels_haemosphere)
col_legend <- list("Cell lineage" = FALSE)
col_anno <- HeatmapAnnotation(
which = "col",
df = col_anno,
col = col_colors,
show_legend = c(FALSE)
)
# Color ramp
color_ramp <- viridis(200, option = "magma")
# Count histogram column annotation
cell_name_count_p <- dplyr::filter(set, tissue == "Progenitor") %>% group_by(fine_labels) %>% dplyr::summarise(n = n()/nrow(.))
cell_name_count_m <- dplyr::filter(set, tissue == "Myeloid") %>% group_by(fine_labels) %>% dplyr::summarise(n = n()/nrow(.))
cell_name_count_p <- dplyr::left_join(data, cell_name_count_p, by = "fine_labels")
cell_name_count_m <- dplyr::left_join(data, cell_name_count_m, by = "fine_labels")
col_bar = HeatmapAnnotation(
"Myeloid fraction" = anno_barplot(cell_name_count_m$n, gp = gpar(border = NA, fill = fill, lty = "blank")),
"Progenitor fraction" = anno_barplot(cell_name_count_p$n, gp=gpar(border = NA, fill = fill, lty = "blank"))
)
hm <- grid.grabExpr(draw(ComplexHeatmap::pheatmap(
mat = m,
name = "Mean delta score",
main = "",
fontsize_row = 12,
scale = "none",
cellwidth = 20,
cellheight = 20,
cluster_rows = FALSE,
cluster_cols = FALSE,
show_rownames = TRUE,
show_colnames = TRUE,
color = color_ramp,
column_split = data$main_labels,
top_annotation = col_bar,
column_title_rot = 90,
border_color = NA)))
return(hm)
}
# Fine label delta score heatmap with fraction of cells histogram annotation
hm_nacl <- hm(subset(so, subset = treatment == "NaCl")@meta.data, color$treatment["NaCl"])
hm_cpg <- hm(subset(so, subset = treatment == "CpG")@meta.data, color$treatment["CpG"])
gridExtra::arrangeGrob(grobs = list(hm_nacl, hm_cpg), nrow = 2, ncol = 1) %>% grid::grid.draw()
ggsave(hm_nacl, filename = "result/plot/progenitor_nacl_haemosphere_deltascore.png", width = 10, height = 5)
ggsave(hm_cpg, filename = "result/plot/progenitor_cpg_haemosphere_deltascore.png", width = 10, height = 5)
Warning message: “Keys should be one or more alphanumeric characters followed by an underscore, setting key from rna_umap_ to rnaumap_” `summarise()` has grouped output by 'fine_labels'. You can override using the `.groups` argument. Warning message: “Keys should be one or more alphanumeric characters followed by an underscore, setting key from rna_umap_ to rnaumap_” `summarise()` has grouped output by 'fine_labels'. You can override using the `.groups` argument.
reduction <- "rna_umap_nno"
options(repr.plot.width = 10, repr.plot.height = 15)
dplot_1 <- DimPlot(so, reduction = reduction, group.by = "rna_snn_res.0.8", label = TRUE) &
ggtitle("Cluster") & xlab("UMAP 1") & ylab("UMAP 2") &
theme(aspect.ratio = 1, legend.position = "none") &
guides(color = guide_legend(ncol = 3, override.aes = list(size = 2)))
dplot_2 <- DimPlot(so, reduction = reduction, group.by = "tissue", label = FALSE) &
ggtitle("Tissue") & xlab("UMAP 1") & ylab("UMAP 2") &
theme(aspect.ratio = 1, legend.position = "bottom") &
scale_color_manual(values = color$tissue, na.value = "dark gray") &
guides(color = guide_legend(ncol = 3, override.aes = list(size = 2)))
dplot_3 <- DimPlot(subset(so, subset = tissue == "Myeloid"), reduction = reduction, group.by = "main_labels", label = FALSE) &
ggtitle("Myeloid") & xlab("UMAP 1") & ylab("UMAP 2") &
theme(aspect.ratio = 1, legend.position = "bottom") &
scale_color_manual(values = color$main_labels_haemosphere, breaks = names(color$main_labels_haemosphere), na.value = "dark gray") &
guides(color = guide_legend(ncol = 3, override.aes = list(size = 2)))
dplot_4 <- DimPlot(subset(so, subset = tissue == "Progenitor"), reduction = reduction, group.by = "main_labels", label = FALSE) &
ggtitle("Progenitor") & xlab("UMAP 1") & ylab("UMAP 2") &
theme(aspect.ratio = 1, legend.position = "bottom") &
scale_color_manual(values = color$main_labels_haemosphere, breaks = names(color$main_labels_haemosphere), na.value = "dark gray") &
guides(color = guide_legend(ncol = 3, override.aes = list(size = 2)))
dplot_5 <- DimPlot(subset(so, subset = tissue == "Myeloid"), reduction = reduction, group.by = "fine_labels", label = FALSE) &
ggtitle("Myeloid") & xlab("UMAP 1") & ylab("UMAP 2") &
theme(aspect.ratio = 1, legend.position = "bottom") &
scale_color_manual(values = color$fine_labels_haemosphere, breaks = names(color$fine_labels_haemosphere), na.value = "dark gray") &
guides(color = guide_legend(ncol = 3, override.aes = list(size = 2)))
dplot_6 <- DimPlot(subset(so, subset = tissue == "Progenitor"), reduction = reduction, group.by = "fine_labels", label = FALSE) &
ggtitle("Progenitor") & xlab("UMAP 1") & ylab("UMAP 2") &
theme(aspect.ratio = 1, legend.position = "bottom") &
scale_color_manual(values = color$fine_labels_haemosphere, breaks = names(color$fine_labels_haemosphere), na.value = "dark gray") &
guides(color = guide_legend(ncol = 3, override.aes = list(size = 2)))
dplot <- dplot_1 + dplot_2 + dplot_3 + dplot_4 + dplot_5 + dplot_6 + plot_layout(ncol = 2)
dplot
ggsave(dplot, filename = "result/plot/dimplot_haemosphere.png", width = 10, height = 15)
Warning message: “Keys should be one or more alphanumeric characters followed by an underscore, setting key from rna_umap_ to rnaumap_” Warning message: “Keys should be one or more alphanumeric characters followed by an underscore, setting key from rna_umap_ to rnaumap_” Warning message: “Keys should be one or more alphanumeric characters followed by an underscore, setting key from rna_umap_ to rnaumap_” Warning message: “Keys should be one or more alphanumeric characters followed by an underscore, setting key from rna_umap_ to rnaumap_”